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FOREWORD 


This  report  was  prepared  by  the  University  of  Dayton  Research  Institute  under  Air  Force 
Contract  No.  F33615-95-D-5029,  Delivery  Order  No.  0004.  The  work  was  administered  under 
the  direction  of  the  Nonmetallic  Materials  Division,  Materials  Directorate,  Wright  Laboratory, 
Air  Force  Materiel  Command,  with  Dr.  James  R.  McCoy  (WL/MLBC)  as  Project  Engineer. 

This  report  was  submitted  in  January  1997  and  covers  work  conducted  from  24  Jan  1996 
through  14  Sep  1996. 


1.  INTRODUCTION 


Development  of  a  damage  progression  analysis  in  composites  contains  an  appropriate 
failure  criterion  and  a  numerical  algorithm  capable  of  incorporating  the  formation  of  new 
boundaries  in  the  material  without  restrictions  imposed  by  initial  mesh  topology.  A  promising 
approach  to  build  such  a  numerical  algorithm  is  local  mesh  overlay  or  local  field  enrichment 
method  described  by  Mote  [1]  and  more  recently  by  Raju  [2],  Fish  [3],  and  Reddy  [4].  The  idea 
is  to  add  additional  degrees  of  freedom  to  the  initial  mesh  via  superimposing  a  patch  of  elements 
at  the  location  of  high  local  stress  gradients,  such  as  those  produced  by  crack  tip,  etc.,  which 
allows  one  to  avoid  changing  the  initial  mesh  and  the  rigidity  matrix.  Instead,  the  new  degrees  of 
freedom  are  added  as  a  separate  block  to  the  solution.  Several  issues  are  not  displayed  in  the 
literature  regarding  this  method;  an  important  question  is  the  accuracy  of  satisfying  the  traction 
continuity  at  the  edge  of  the  patch. 
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2.  PROBLEM  FORMULATION 


Consider  an  elastic  body  occupying  a  volume  V,  containing  a  crack  S  (Figure  la). 
Displacement  and  traction  boundary  conditions  are  imposed  over  surfaces  dVu  and  dVr, 
respectively,  where  dVu  +  dVr,  =  dV.  We  shall  seek  the  displacements  as  a  superposition  of  two 
terms: 


Ui  (x,y,z)  =  uf{x,y,z)  +  uf(x,y,z),  ( 1 ) 

where  uf(x,y,z)  are  functions  continuous  through  the  entire  body  V,  and  u^{x,y,z)  are  functions 
discontinuous  at  the  surface  S.  Displacement  field  (1)  is  assumed  to  be  kinematically  admissible. 
In  addition  we  define  a  volume  F,  which  includes  the  crack  surface  S,  which  also  may  be  a  part 
of  the  boundary  9r.  The  functions  uf(x,y,z)  are  required  to  vanish  outside  and  at  the  boundary 
of  the  volume  F 


uf{x,y,z)  =  0,  [x,y,zl {.x,y,z)  €  F  and  {x,y,z)  e  ^-5}.  (2) 

It  is  noted  that  no  assumptions  are  made  by  introducing  Equations  (1)  and  (2).  The  stress  fields 
corresponding  to  the  displacement  fields  wf  (x,y,z)  and  in  general  will  be  discontinuous  at  the 
boundary  9r  and 

lim +  pn)nj  =  lim  [<T”(M - pn)  +  oUM - pn)]nj,  P>0  (3) 

)S->0  ■'  p-^0 

where  M  e  3F  and  n  is  an  outside  normal  to  3F.  It  is  of  practical  interest  of  how  to  obtain  the 
class  of  solutions  when  (7y{M)  =  0,  Af  e^F .  The  importance  of  this  requirement  can  be 
appreciated  by  using  numerical  approximations  to  obtain  the  displacement  functions  in  Equation 
(1).  If  this  requirement  is  not  satisfied,  then  the  stress  <jy  will  experience  a  discontinuity  at  the 
boundary  9F  which,  if  not  properly  built  into  displacement  approximation  uf(x,y,z),  will  cause 
oscillatory  behavior.  Yet  the  attractiveness  of  the  entire  concept  of  patch  superposition  lays  in 
the  ability  to  use  the  same  approximation  for  uf(x,y,z)  as  for  the  uncracked  body.  The 
following  variational  equation  will  be  shown  to  provide  the  condition  <Jy(M)  =  0,Medr: 


2 


(a) 


Figure  1 .  Illustration  of  Displacement  Field  Superposition  Concept:  Multilayered  Plate 
with  (a)  a  Hole  and  (b)  a  Crack. 
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-IJi  v(4.n + 4j)'xi'' 

V  v-r  r 

+jj<7lnjufds 


dr 


(4) 


JJ  Tf^ufds  +  Jj  [T^iuf^  +uf)-  Triuf-  +  uf)]ds 


dVr 


=  0 


where  =  ^(«,j  +  CTy  =  Qijki’^(k,iy  ^(“(«J))  =  ^Qijki^(ij)^ik,i)’  Qijki  are  the  elastic 

constants,  superscripts  “c+”  and  “c-”  denote  the  two  crack  faces,  and  T^aad  Tf  denote  the 
external  tractions  applied  to  the  latter  surfaces.  The  variation  upon  uf{x,y,z)  yields: 


111  4i.i^Uv + Ijj  Kj + + 

V-r  r 

+Jj  ((er^l  v-r  ~  (y^lrYjSufds  + 

dr 

JJ (afjnj  -  Tf^)Su^,ds+ll (af  -  alpnj  -T*  +  Tr]5uf  =  0 
dVr  s 


(5) 


The  surface  integral  upon  9r  is  obtained  by  taking  into  account  the  respective  surface  integral  in 
(5),  where  n,-  is  an  outwards  relative  to  volume  F  normal  vector.  The  traction  contributions  of  the 
stresses,  created  by  the  uf  (x,  y,z)  displacement  field  coming  from  the  first  and  second  volume 
integral  in  Equation  (4),  are  distinguished  because  according  to  Equation  (3),  they  can  be  unequal 
on  the  boundary  9r.  However  due  to  arbitrary  variation,  5uf  on  the  boundary  dF,  we  have 

((cr9.|^_r(M)-c7g.|r(M))n,-  =0,Me5F  (6) 


The  surface  integrals  upon  the  crack  surface  and  the  external  surface  dVr  are  responsible  for 
satisfaction  of  the  applied  external  traction  boundary  conditions. 
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Variation  upon  Mf(jc,}',z)  provides: 

r 

+ jj  (cr°|  r  +  G^i^Hjduids  +  Jj  Sa^jrijufds 

dr  dr 

+JJ ((erf  +  al)nj  -  +  JJ ((trf  +  c^^rij  -  Tr)5urds  =  0 

s  s 

Both  surface  integrals  over  3r  vanish,  one  due  to  the  boundary  condition  (2)  and  the  second  due 
to  G‘-j(M)nj  =  0,M  e  dr ,  which  follows  from  Equation  (6). 
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3.  COMPOSITE  LAMINATE  WITH  A  HOLE  CONTAINING  MATRIX  CRACKS 

Consider  a  rectangular  orthotropic  plate  containing  a  circular  hole  having  a  diameter  D  as 
shown  in  Figure  lb.  The  plate  consists  of  N  plies  of  total  thickness  H  in  the  z-direction  and  has  a 
length  L  in  the  x-direction  and  width  A  in  the  y-direction.  Uniaxial  loading  of  the  plate  in  the 
x-direction  is  considered.  At  the  opposite  edges  of  the  plate,  x=0,L,  constant  displacement  in  the 
x-direction  is  prescribed,  and  other  displacement  components  at  these  edges  are  presumed  to  be 
zero; 

u^(0,y,z)  =  -ui,  M  (0,y,z)  =  uJ0,y,z)  =  0 

(7) 

u^(L,y,z)  =  ui,  UyiL,y,z)  =  u^{L,y,z)  =  0. 


A  cylindrical  coordinate  system  r,9,z  is  introduced  at  the  center  of  the  hole.  The  0=0°  coincides 
with  the  x-axis. 


A  crack  of  length  /  emanating  from  the  hole  edge  at  6=Yo  in  the  s-th  ply  and  propagating 
in  the  direction  0=a  is  considered.  The  crack  surface  S  is  defined  as 


X  =  -cosy/Q+^-lcosa  +  Xc;,  0  <  |  <  1, 

3 

y  =  — sin^o  -H^-Zsina  +  y^,  0  <  |  <  1, 


(8) 


where  the  s-th  ply  occupies  a  region  z^^  <  z  < 
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4.  SPLINE  APPROXIMATION  OF  DISPLACEMENTS 


The  displacement  field  uf(x,y,z)  is  approximated,  according  to  larve  [5].  Cubic  spline 
functions  are  used.  The  displacements  are  continuous  through  the  laminate;  the  strains  and 
stresses  are  continuous  within  a  homogeneous  ply.  Curvilinear  transformation,  which  maps  the 
x,y  plane  of  the  plate  with  a  hole  into  a  region  0<p<l,  0<<j)<27t,  was  defined  as  follows: 

X  =  — Fj(p)cos0  +  L-F2(p)a(^)  +  x^ 

I  .  (9) 

y  =  —Fiip)sm(j)  +  A-  F2  (p)y3(^)  + 


where 


Fi(p)  = 


\l  +  K-p,p<Ph 
(l  +  K-pf^)(l-p) 


1-P/, 


,  P;,  <  P  <  1 


F2(p)  = 


{^’p^ph 

P-Ph 


^~Ph 


,  Ph<P<l 


The  coordinate  line  p=0  describes  the  contour  of  the  hole,  and  the  coordinate  line  p=l  describes 

the  rectangular  contour  of  the  plate.  Inside  the  near-hole  region,  D/2<p<(l+K)D/2,  which 

corresponds  to  0<p<Ph,  a  simple  relationship  between  the  cylindrical  coordinates  and  the 

D  Dk 

curvilinear  coordinates  p,(l)  exists:  p - = - p  and  0  =  0.  The  width  of  this  region  will  be 

2  2 


chosen  three-hole  radii,  i.e.,  Kph=3.  Beyond  this  region  a  transition  between  the  circular  contour 
of  the  opening  and  the  rectangular  contour  of  the  plate  occurs.  Functions  a((j))  and  ^(0)  were 
defined  [5],  so  that  the  parametric  equations  x=a(0)+Xc,  y=P(0)+yc  describe  the  rectangular 
contour  of  the  plate,  where  corresponds  to  0<x^,  y=A,  corresponds  to  x=0, 

0<y^,  corresponds  to  0<x<L  and  y=A  and  0^^^<0<27i  corresponds  to  x=L,  0<y<A. 


Displacement  components  in  the  x,  y  and  z  coordinate  directions  were  approximated  by 
using  cubic  polynomial  spline  function  of  generalized  curvilinear  coordinates  p,  0,  and  z  as 
follows: 


H, = (10) 
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where  unknown  vectors  Ug,  V^,  Wj  contain  unknown  displacement  spline  approximation 
coefficients.  Bold  type  here  and  below  will  be  used  to  designate  matrices  and  vectors,  and  the 
superscript  star  means  the  transpose  operation.  The  vector  of  the  three-dimensional  spline 
approximation  functions  was  defined  as: 

\  ^  ^  3 

q  =  l  +  (j-l)(nj+3)  +  ii-V)(nj.  +  3)k,l  =  l,...,n^  +  3J  =  =  l,...,m  +  3. 


where  sets  of  B-type  cubic  basis  spline  functions  along 

each  coordinate  were  built  upon  subdivisions:  0  =  po  <  Pi  < ...  <pm  =  1, 0  =  <  —  <  ^*  = 

271:,  =  zo<Zi< ...  <  z„^  =  z^^^  so  that  the  s-th  ply  occupies  a  region  <z^  z^^\  and  n,  is  the 

number  of  sublayers  in  each  ply.  The  subdivision  of  the  p  coordinate  is  essentially  nonuniform. 
The  interval  size  increases  in  geometric  progression  beginning  at  the  hole  edge.  The  region  0  <  p 
<  Ph  in  which  the  curvilinear  transformation  is  quasi-cylindrical,  is  subdivided  into  mo  intervals, 
so  that  pf^  =  p„^ .  Numbers  of  intervals  of  subdivision  m,  k,  ns  in  each  direction,  along  with  the 
mesh  nonuniformity  characteristics  such  as  mo  and  the  consecutive  interval  ratio,  determine  the 
accuracy  of  the  solution  and  the  size  of  the  problem.  Nonsquare  boundary  matrices 

and  constant  vectors  Eo,El  are  defined,  so  that  the  approximation  (10)  provides  a 
kinematically  admissible  displacement  field  for  any  coefficients  Us,Vs,Ws,  i.e.,  satisfying 
boundary  conditions  (7).  The  components  of  vectors  Eo,El  are  equal  to  1  if  the  same 
components  of  nonzero  at  p=l,  (|)^*^<(j)<(|)^^^  (x=0, 0<y<A)  and  p=l,  (|)^^^<(l)<27t  (x=L, 

0<y<A),  respectively.  All  other  components  of  the  vectors  Eo,El  are  equal  to  zero.  The 
boundary  matrices  are  obtained  by  deleting  a  number  of  rows  from  the  unit  matrix.  The  rows  10 
deleted  are  the  ones  having  a  nonzero  scalar  product  with  Eq  or  El. 


4.1  Model  1:  Spline  Approximation  of  Displacements  in  the  Overlay  Patch 

The  overlay  displacement  field  ufix,y,z)  was  also  approximated  using  cubic  spline 
functions.  The  location  of  the  overlay  mesh  was  defined  by  equations: 


D 


X  =  — cosy/  + 1  •  /  +  — (cosCv^o  +  «)  -  cos(y/  -  a)) 

2  \  2  j 


cosa  +  x.. 


y  =  —smy/  +  ^-(l+— (cos(  i/^’o  +  a)-  cos(y/  -  a))' 
2  V  2  J 


sma-Hy^, 


(12) 
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where  0<^<1,  \i/o-Ai<\i/<\|/o+A2.  The  boundary  ^=0  coincides  with  the  part  of  the  hole  boundary, 
and  the  boundary  ^=1  is  a  straight  line  going  through  the  crack  tip.  Boundaries  \)/=  xj/q-Ai,  \|/= 
\};o+A2  are  parallel  to  the  crack  face.  Subdivisions  are  introduced  0=^o<^i<^2<  —  =1’ 

i/Tq  -  A|  =  <  Y-ki+i  ¥o  <  ¥i  <•••<  Wk2  -¥o'^  ^2-  subdivisions  are  uniform, 

and  the  crack  surface  \|/=\|/o  is  a  coordinate  line  of  the  subdivision,  total  number  of  intervals  of 
the  \|/  subdivision  is  kc=ki+k2.  The  coordinate  z  is  subdivided  into  the  same  sublayers  as  the 
laminate,  i.e.,  =  Zo  <  <•••<  example  of  the  initial  mesh  (9)  with  m=18, 

mo=14,  q=1.0,  Kph=3,  k=48,  L/A=2,  L/D=10,  Xc=L/2,  yc=A/2  is  given  in  Figure  2a,  an  overlay 
mesh  with  (x=90°,  v|/o=90°,  mc=9,  ki=k2=7,  Ai=A2=30°,  1=0. IL  is  shown  in  Figure  2b  and 
superimposed  meshes  in  Figure  2c.  Basic  systems  of  spline  functions  are  built  upon  the  above 
subdivisions.  Two  independent  sets  of  splines  are  built  upon  \j/:ki+3  functions  on  the  interval 
y/Q-Ai  =  <  V'-jtj+i  <...<  V^o  ^2+3  functions  on  the  interval 

¥o  <  ¥i  <■■■<  Wk2  =  V^o  +  ^2-  "That  way  no  continuity  conditions  are  imposed  at  \i/=\|to-  Spline 
approximation  of  the  overlay  displacements  can  be  expressed  as 

ui = c^^cu:,  = cixX>  < = oixc^:.  (13) 

{Xc}q  = 

q  =  /  +  (7-l)(n5  +  3)  +  (i-l)(n^  +  3)[^^  +  6],  I  =  +  =  \,...,k^  +  6,i  =  l,...,m^  +3. 

Vectors  Uc,Vc,Wc  contain  the  unknown  spline  approximation  coefficients.  Sets  of  basic  spline 
functions  are  denoted  as  (‘^)  j  » I*®*/  I z-functions  are  the  same  as 

in  (10).  To  provide  the  boundary  conditions  (2),  the  boundary  matrices  have  to  be  defined,  so 
that  all  components  of  the  spline  function  vector  different  from  zero  at  \}t=\|/o-Ai,  V|/=\}ro+A2,  ^=1 
and  z=z^®‘’\  z=z^*^  will  not  contribute  to  displacement  values.  According  to  boundary  properties 
of  spline  functions  given  by  larve  [5],  the  only  components  of  the  spline  function  vector 
contributing  at  these  surfaces  are  the  ones  containing  functions  j  =  1,  k(,  +  6; 

i  =  mj.  +  3  and  Z/(z),  /  =  1,  +  3 .  The  boundary  matrices  are  obtained  from  unit 

matrices  by  deleting  the  rows  having  a  nonzero  scalar  product  with  E‘^,  where  a  component  of 
is  equal  to  unity  if  the  corresponding  component  of  the  vector  Xc  contains  the  aforementioned 

functions  and  zero  otherwise. 
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Figure  2.  (a)  Original  Mesh,  (b)  Overlay  Mesh,  and  (c)  the  Superposition. 
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4.2  Model  2:  Spline  Basis  Decomposition  Method 

An  alternative  approach,  which  circumvents  the  compatibility  conditions  resolved  by 
functional  (4),  is  to  create  a  complete  system  of  basic  functions  in  the  spline  function  space, 
which  provides  displacement  discontinuity  along  the  crack  surface  S  described  by  Equation  (8). 

It  is  possible  to  build  such  a  basis  without  modifying  the  spline  functions  of  approximation  (10) 
but  rather  by  adding  locally  new  splines  to  describe  the  discontinuity  of  displacements.  The  idea 
will  be  illustrated  on  a  one-dimensional  example.  Consider  an  interval  subdivided  into  seven 
subintervals  (Figure  3).  We  shall  build  the  spline  approximation  of  a  function  f(x),  which  is 
discontinuous  at  the  point  x=X4,  and  continuously  differentiable  elsewhere  on  the  interval  [xo,X7]. 
A  straightforward  solution  lies  in  using  two  sets  of  spline  basis  functions  built  over  intervals 
[xo,X4]  and  [X4,X7],  respectively.  The  total  number  of  spline  functions  will  be  m+6,  where  m  is 
the  number  of  intervals  (m=7).  In  this  case,  however,  most  of  the  splines  of  such  a  basis  do  not 
coincide  with  the  spline  functions  of  continuously  differentiable  approximation  on  the  entire 

interval  [xo,X7]  as  shown  in  Figure  3a.  These  functions  are  denoted  {X3i(x)}i=i  „,+3,  so  that  a 
function  f 


m+3 

/(x)= 


1=1 


is  continuously  differentiable  at  each  point  between  Xq  and  X7.  New  spline  functions  X’3j(x)  are 
formed  by  partitioning  splines  which  X3  i(x4)>0.  These  new  splines  are  defined  as  follows 


X3,(A;),a:>A:4 

0,X<X4 


In  the  present  example  there  are  three  new  splines  created  i=5,6,7  as  shown  in  Figure  3b.  A 
function  f(x) 

m+3  7 

f(x)  =  5;  +  X  Ji'X3,iU)  (14) 

i=l  i=5 

is  discontinuous  at  x=X4  and  continuously  differentiable  at  all  other  points  of  the  interval.  It  can 
also  be  shown  that  functions  {X3,i(x)}i= 1,^+3  and  {X’34(x)}i=5j  form  a  complete  set  of  basis 
functions.  Thus  we  have  built  a  complete  set  of  basis  functions  by  adding  several  new  splines  to 
the  initial  continuous  approximation  spline  functions. 
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Figure  3.  Spline  Basis  Decomposition  Method  on  a  One-Dimensional  Example:  (a)  Original, 
Continuously  Differentiable  Basis  Functions  and  (b)  Additional  Spline  Functions. 


In  the  case  of  a  crack  emanating  from  the  hole  edge,  we  need  to  build  essentially 
a  two-dimensional  set  of  splines  analogous  to  the  X’3  i(x)  in  Equation  (14).  In  the  z-direction  the 
same  spline  functions  as  in  the  previous  section  will  be  used.  A  parametric  representation  of 
the  x,y  plane  is  defined: 


X  =  — cosy/Q  +  |-/cosa-Tsina  +  Xp,(|  >  0 

'2d 


y  =  — siny^Q  -l-^•/sina-Tcosa  +  y^,^  >  0 
2 


The  spline  functions  in  Equation  (1 1)  are  examined,  and  those  nonzero  at  the  crack  line  t=0  are 
partitioned  to  create  new  spline  functions 

{x;}^=j;;(p)4./«z!‘>a),  (16) 


fi;(p)<i>/«z,W(z)= 


J  Riip)^jWZl^\z),r>0 

1  0,T<0 
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Additional  or  superimposed  spline  approximations  can  be  written  similar  to  Equation  (13): 


K  =  CSA. =  cixX. =  cixM- 


where  the  vector  of  spline  functions  is  defined  similar  to  (1 1),  but  only  the  spline  functions  (16) 
are  used.  Variational  formulation  in  this  case  is  straightforward,  since  no  incompatibility  exists 
between  shape  functions  (11)  and  (17).  The  variational  equation  will  be  simplified  as: 


V  V 

jj  +«?)-  7r(»r +“?)] 
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5.  NUMERICAL  RESULTS 


A  unidirectional  AS4/3501-6  laminate  [9O2]  was  considered  first.  The  mechanical 
properties  were  Ei  =  138  GPa,  E2  =  E3  =  10.3  GPa,  V13  =  V12  =  0.3,  V23  =  0.55,  G12  =  G13  =  5.52 
GPa  and  G23  =  3.45  GPa.  The  geometric  dimensions  were  L  =  6.35  cm,  A  =  L/2,  and  the 
diameter  of  a  central  hole  was  D  =  L/10.  A  crack  /  =  0.1  L  was  considered  emanating  from  the 
hole  edge  at  0  =  90°  (\|/o  =  90°)  in  the  a  =  90°  direction.  The  p,(|)  subdivision  with  m  =  18,  mo  = 
14,  q  =  1.0,  Kph  =  3,  and  k  =  48,  shown  in  Figure  2a,  along  with  two  sublayers  in  the  z-direction, 
were  utilized  to  approximate  the  u  ®(x,y,z)  displacement  field. 

Model  1  was  used  to  approximate  the  overlay  displacements.  The  overlay  mesh  was 
imposed  only  on  one  side  of  the  crack,  so  that  Ai  =  30°  and  A2  =  0,  with  10  \|/-intervals  between 
\l/o-Ai  and  \|/o.  Ten  uniform  ^  intervals  were  used.  Stress  distribution  in  the  cross  section  9  =  90° 
as  a  function  of  distance  from  the  hole  edge  is  examined  in  Figure  4a.  The  Gxx  stress  component 
in  the  plate  without  a  crack  is  shown  for  comparison  by  a  thick  dashed  line.  The  stress  results  are 
normalized  to  the  average  far-field  stress  in  an  uncracked  laminate  in  the  x-direction  ,  calculated 
as 

HA 

^0  =  \\<y  xx(<Uy,z)dydz. 

00 

The  results  for  the  cracked  laminate  are  obtained  using  the  variational  equation  (4),  neglecting 
the  surface  integral  upon  the  boundary  dP.  Stress  distribution  in  the  cracked  laminate  is  shown 
by  the  solid  and  thin  dashed  lines.  In  addition  the  crack  opening  displacement  is  also  shown, 
normalized  to  the  applied  displacement  Ul  [Equation  (7)].  Along  the  crack  surface  in  the  vicinity 
of  the  crack  tip  CTxx,  stress  exhibits  oscillatory  behavior.  However  the  true  effect  of  the 
incompatibility  of  the  overlay  and  underlying  meshes  is  clearly  illustrated  in  Figure  4b,  where  the 
stresses  along  the  boundary  of  the  patch  region  \|f  =  \|/o-Ai  are  shown.  The  thin  dashed  line  shows 
the  stress  designated  in  (5)  as  cr^|  ,  and  the  solid  line  stress  cr^|r(M)  +  where 

Me  3r.  A  significant  stress  discontinuity  is  obvious.  The  solid  dashed  line  shows  the  Cxx  stress 
in  a  cross  section  symmetric  to  v|/  =  Xj/o-Ai  against  the  plane  6  =  90°.  Ideally,  all  three  curves 
shown  in  Figure  4b  have  to  coincide. 

Model  2  will  be  used  to  approximate  the  overlay  displacements  in  the  same  problem.  The 
displacements  u  ®(x,y,z)  are  the  same  as  before.  The  superimposed  shape  functions  are  formed 
according  to  rule  (16)  and  no  additional  mesh  is  required.  Variational  Equation  (18)  is  used. 
Figure  5a  shows  the  CJxx  stress  distribution  in  the  cross  section  6  =  90°.  The  thick  dashed  line 
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(a) 


(b) 


Figure  4.  Mesh  Overlay  Model  1,  at  (a)  0  =  90°  and  (b)  Overlay  Mesh  Boundary. 
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(a) 


(b) 


Figure  5.  Mesh  Overlay  Model  2,  at  (a)  0  =  90°  and  (b)  at  Model  1  Overlay  Mesh  Boundary. 
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was  copied  from  Figure  4a.  The  stress  calculated  by  using  Model  2  for  overlay  displacement 
approximation  is  shown  by  a  thin  solid  line.  No  oscillations  of  the  cTxx  stress  along  the  crack  face 
is  present.  There  is  a  region  in  the  vicinity  of  the  crack  tip  marked  as  a  crack  tip  boundary 
condition  error.  In  this  region  both  the  crack  opening  displacement  and  the  stress  normal  to  the 
crack  faces  are  nonzero.  By  increasing  the  density  of  the  subdivision  by  splitting  four  intervals 
near  the  crack  tip,  this  region  reduces.  The  result  for  the  increased  subdivision  is  shown  by  a 
thick  solid  line.  The  stresses  at  the  same  cross  sections,  as  shown  in  Figure  4b,  are  displayed  in 
Figure  5b.  The  stress  obtained  by  using  Model  2  along  the  cross  section  \|r  =  \|/o-Ai,  and  the  one 
symmetric  against  the  plane  6  =  90°,  are  indistinguishably  close-solid  lines.  The  dashed  line  has 
been  copied  from  Figure  4b  and  shows  overall  agreement  with  the  Model  2  results,  except  a  kink 
near  the  end  of  the  pitch  which  is  likely  the  result  of  untreated  compatibility  of  the  meshes  in 
Model  1.  Figure  6  shows  the  hoop  and  the  radial  stress  in  uncracked  and  cracked  laminates 
obtained  using  Model  2.  The  displacement  field  in  the  uncracked  laminate  was  approximated 
according  to  (1 1)  only.  Stress  is  plotted  versus  the  0  angle  at  the  midsurface.  The  presence  of 
the  crack  at  0  =  90°  increases  the  hoop  stress  at  0  =  270°.  An  important  observation  is  that  the 
accuracy  of  satisfaction  of  the  traction-free  boundary  condition  at  the  hole  edge  does  not 
deteriorate  while  enriching  the  spline  approximation  basis  with  functions  (16).  Indeed,  the 
stress  distribution  for  uncracked  and  cracked  laminates  is  practically  the  same. 


Figure  6.  Radial  and  Hoop  Stresses  at  the  Midsurface  Around  the  Circumference  of  the  Hole  of 
the  Uncracked  and  Cracked  Laminate. 
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A  [O2]  laminate  is  considered  next.  A  matrix  crack,  tangential  to  the  hole  edge,  is 
modeled  as  two  cracks  of  length  /  =  0.1  L  emanating  at  0  =  90°  (\|/o  =  90°)  in  a  =  0°  and  a  =  180° 
directions.  All  of  the  following  results  are  obtained  with  Model  2  overlay  displacement 
approximation.  The  applied  displacement  was  Ul  =  Uq  =  6.95  x  10'^  m.  Contour  plot  of  the  Uy  in 
units  of  meters  is  shown  in  Figure  7b  inside  the  region  designated  in  Figure  7a.  The  Uy  is 
positive  above  the  crack  and  negative  below  it,  reflecting  the  fact  that  this  crack  is  open.  The 
normal  Cyy  and  shear  stresses  are  shown  in  Figure  7c  and  7d.  The  presence  of  the  crack 
relieves  both  stresses  in  the  area  above  the  crack  with  only  stress  concentrations  at  the  crack  tips. 

Finally,  a  [0/90]s  laminate  is  considered.  The  same  geometric  dimensions  as  before  are 
used.  Each  ply  is  subdivided  into  two  sublayers.  Two  matrix  cracks  \|/o  =  90°  and  \|/o  =  270°  are 
considered  in  the  90°  ply,  as  shown  in  Figure  8a.  The  o^x  stress  around  the  hole  edge  is 
examined  at  different  through-the-thickness  locations.  Figure  8b  displays  the  stresses  at  the 
midsurface  and  the  0/90  interface  in  the  90°  ply.  The  dashed  lines  are  used  for  uncracked 
laminates,  and  the  solid  lines  for  cracked  laminates.  At  0  =  90°  and  270°  locations  at  the 
midsurface,  Cxx  vanishes  at  the  crack  surface.  This  stress  redistribution  is  very  local,  and  at 
approximately  ±15°  away  from  the  crack  locations,  the  stress  is  equal  to  that  in  the  uncracked 
laminate.  At  the  ply  interface  the  crack  produces  higher  stresses  because  it  is  impending  the 
interface.  In  the  0°  ply  the  stresses  at  0  =  90°  and  270°  are  slightly  higher  in  the  cracked  laminate 
than  the  uncracked  one.  It  is  expected  due  to  the  loss  of  load  carrying  capacity  in  the  90°  ply. 
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The  window  in 
which  the  contour 
plots  are  performed 


Polar  Angie 


cracks  in  90°  ply 


(c) 


Figure  8.  Crack  Location  and  the  (a)  Coordinate  Systems;  Hoop  Stress  Redistribution  in  the 
(b)  90°  Ply  and  (c)  0°  Ply  due  to  Matrix  Cracking.  Solid  lines  =  cracked  specimen 
dashed  lines  =  uncracked  specimen. 
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